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Recent progress 

We have previously reported work related to basic technique development in phase 
unwrapping and generation of digital elevation models. In the final year of this work we 
have applied our technique work to the improvement of DEM’ s produced by SRTM. In 
particular, we have developed a rigorous mathematical algorithm and means to fill in 
missing data over rough terrain from other data sets. 

We illustrate this method by using a higher resolution, but globally less accurate, DEM 
produced by the TOPSAR airborne instrument over the Galapagos Islands to augment the 
SRTM data set in this area. We combine this data set with SRTM to use each set to fill in 
holes left over by the other imaging system. The infilling is done by first interpolating 
each data set using a prediction error filter that reproduces the same statistical 
characterization as exhibited by the entire data set within the interpolated region. After 
this procedure is implemented on each data set, the two are combined on a point by point 
basis with weights that reflect the accuracy of each data point in its original image. 

In areas that are better covered by SRTM, TOPSAR data are weighted down but still 
retain TOPSAR statistics. The reverse is true for regions better covered by TOPSAR. 

The resulting DEM passes statistical tests and appears quite feasible to the eye, but as this 
DEM is the best available for the region we cannot fully verify its accuracy. Spot checks 
with GPS points show that locally the technique results in a more comprehensive and 
accurate map than either data set alone. 

The results are described in a paper submitted to IGARSS-2004, which is attached to give 
a complete description of the method and its results. The abstract for this work follows: 

“We present a method for correcting large-scale artifacts and for interpolating regions of 
missing data in TOPSAR DEMs. Artifacts are eliminated by subtracting the difference 
between the TOPSAR DEM and the SRTM DEM. The inversion method using a 
prediction error (PE) filter was used to interpolate the missing data holes so that the 


interpolated regions would have the same spectral content as the originally valid regions 
of the TOPSAR DEM. Along with the PE filter, the SRTM DEM for the same area is 
used as another constraint in the interpolation. We used cross-validation to obtain the 
optimal weighting of the PE filter and SRTM DEM constraints.” 

In previous years, we have been able, using network flow techniques, to unwrap very 
large (100 km x 800 km) ERS interferograms , which were illustrated in large scale maps 
of Alaska. We also modified the statistical description of the network flow paths to allow 
the methods to be used for deformation interferograms, and applied this version to 
unwrapping the Hector Mine interferogram. The method has proved exceptionally stable 
and accurate, and is described in detail in a Ph. D. thesis by Curtis Chen. 

Results presented in several journal articles as listed below and in Chen’s thesis, as listed 
in his abstract: 

“Two-dimensional phase unwrapping is a key step, and often the most significant error 
source, in the analysis of synthetic-aperture-radar interferograms. In the interferometric 
technique, very accurate measurements of the Earth’s topography or its surface 
deformation are derived from radar-image phase data. Phase, however, is defined only 
modulo 2 pi rad, so a resulting 2-D array of measurements is wrapped with respect to 
some modulus or ambiguity. These data must be unwrapped to provide meaningful 
information. For this purpose, we introduce a new, nonlinear constrained-optimization 
approach in which i) defined cost functions map particular unwrapped solutions to scalar 
costs and ii) a solver routine computes minimum-cost solutions. Previous efforts have 
focused mainly on simple cost functions that have yielded efficient — but not necessarily 
accurate — algorithms. These inaccuracies seriously degrade the effectiveness of the 
interferometric technique and can preclude useful geophysical interpretation of the data. 
We propose a new set of nonconvex, statistically based cost functions through which we 
treat phase unwrapping as a maxim um a posteriori probability estimation problem. That 
is, we derive approximate, application-specific statistical models for the problem 
variables. Based on these models, we cast phase unwrapping as an optimization problem 
whose objective is to find the most physically probable unwrapped solution given the 
observable quantities: wrapped phase, image intensity, and interferogram coherence. We 
prove that the resulting problem is NP-hard, and we develop nonlinear network-flow 
solver techniques for approximating solutions to this problem. Extending our statistical 
framework and network methods, we also present a tiling heuristic for applying our 
algorithm to large data sets. Performance tests on topographic and deformation data 
acquired by the ERS-1 and ERS-2 satellites suggest that our algorithm yields superior 
accuracy and competitive efficiency as compared to other existing algorithms.” 

We have in addition published three journal articles and several conference presentations 
describing various aspects of this problem. 

This study has produced one Ph.D. thesis (Curtis Chen, now at JPL) and will be the basis 
of several chapters in Sang-Ho Yun’s Ph.D. dissertation. 
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Abstract — We present a method for correcting large-scale 
artifacts and for interpolating regions of missing data in TOPSAR 
OEMs. Artifacts are eliminated by subtracting the difference 
between the TOPSAR DEM and the SRTM DEM. The inversion 
method using a prediction error (PE) filter was used to interpolate 
the missing data holes so that the interpolated regions would have 
the same spectral content as the originally valid regions of the 
TOPSAR DEM Along with the PE filter, the SRTM DEM for 
the same area is used as another constraint in the interpolation. 
We used cross-validation to obtain the optimal weighting of the 
PE filter and SRTM DEM constraints. 


I. Introduction 

The recent SRTM mission produced worldwide topographic 
data at 90m postings. For many analyses, however, finer- 
scale elevation data are required. We are currently modeling 
volcanic processes in the Galapagos Islands, and in fact have a 
lOm-posting DEM of Isabela and Femandina acquired by the 
TOPSAR airborne instrument. TOPSAR DEMs are produced 
from cross-track interferometric data acquired with NASA’s 
AIRSAR system mounted on a DC-8 aircraft Although the 
TOPSAR DEMs have a higher resolution than other existing 
data, they sometimes suffer from missing data holes and other 
artifacts due to layover (Fig. 2a), flight planning limitations, 
and roll of the aircraft (Fig. 1). Fortunately, die SRTM DEM 
has fewer missing data than the TOPSAR DEM, and thus the 
former provides some information on the missing regions of 
the latter. 

TOPSAR and SRTM missions are the primary sources 
for DEMs derived from single-pass interferometric data. 
However, differences in their system parameters such as 
altitude and swath width (Table I) have resulted in slightly 
different products. In particular, TOPSAR DEMs have better 
resolution, while SRTM DEMs have better data coverage and 
consistency. In this presentation, we describe a method for 
combining the two DEMs to produce a new DEM that has 
a resolution of TOPSAR DEM and a consistency of SRTM 
DEM 


TABLE I 

TOPSAR MISSION VS. SRTM MISSION 


Mission 

TOPSAR 

SRTM 

Platform 
Nominal altitude 
Swath width 
Baseline 
DEM resolution 
DEM coord, system 

DC-8 aircraft 

9 km 

10 km 
2.583 m 

10 m 
none 

Space shuttle 
233 km 
225 km 
60 m 
90m 
Lat/Lon 


II. Method 

A. TOPSAR DEM Registration 

The original TOPSAR DEM was not registered to a geo- 
referenced frame. Thus, we registered the TOPSAR DEM to 
the SRTM DEM, which was registered in a latitude/longitude 
coordinate system. This registration was done with one affine 
matrix and one translation vector. We picked 10 tie points from 
each DEM and solved for the six unknowns - four from the 
affine matrix and two from the translation vector in a least- 
squares sense. Then we resampled the TOPSAR DEM such 
that the resulting DEM would have a pixel spacing exactly 
nine times smaller than the pixel spacing of the SRTM DEM. 
Cross-correlation was used to check the remaining offset. 

B. Artifact Elimination 

An artifact most likely due to roll of the DC-8 aircraft 
was found in the TOPSAR DEM (Fig. 1). We averaged the 
TOPSAR DEM by taking 9 looks in both east and north 
directions and subtract it from the SRTM DEM. With the 
maximum amplitude of over 20 meters, this artifact would 
result in substantial errors in many analyses. For example, 
if this DEM is used for DInSAR application with an ERS2 
data pair that has a perpendicular baseline of about 400 
meters, the resulting interferogram would contain one fringe of 
spurious deformation. The space shuttle data does not exhibit 
significant roll artifacts as its swath is much larger than that 
of TOPSAR (Table I). Thus, the roll artifact can be eliminated 
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Fig. 2. Example images of (a) TOPSAR DEM, and (b) SRTM DEM of the 
same area. The grayscale is in meters 


Fig. 1. Artifact due to rolling of the aircraft during the TOPSAR mission. 
The grayscale is in meters 


by subtracting the difference between the TOPSAR DEM and 
the SRTM DEM. 

C. Prediction Error Filter 

Based on an assumption that the regions missing data have 
the same spectral content as the regions with valid data, we 
generate a PE filter [1] such that it would nullify the valid 
regions of the TOPSAR DEM. Given this PE filter, we solved 
for die missing regions such that the interpolated region would 
also be nullified by the PE filter. 

The PE filter, Fp#, is generated such that it minimizes the 
following objective function, 

I |Fp£ Xexiating ~ 0| | 0) 

where x existing is the valid data from the TOPSAR DEM, 
and 0 is a null vector that has the size of Existing • These 
matrix and vector expressions are used just to indicate their 
linear relationship. All inverse problems in this study were 
solved with conjugate gradient method, where forward and 
adjoint functional operators are used instead of the explicit 
form of matrices. Otherwise M x N times larger memory 
would be required to store the matrix, where M x N is the 
size of the image. The output of each operator becomes an 
input parameter of a conjugate gradient solver at each iteration. 

D. Interpolation 

Originally, the SRTM mission produced topographic data 
at 30m postings (1 arc second). However, those data sets 
outside of the United States are not available to the public yet. 
Intead, DEMs at 90m postings (3 arc second) are available to 
download. These data were produced by taking 3 looks in both 
east and north directions. In a similar way, we start with the 
assumption that the pixel value in an SRTM DEM with 90m 


postings is equivalent to the average value of a 9 by 9 pixel 
window centered at the same location in the TOPSAR DEM. 
We introduce this averaging process in our inverse problem 
with the hope that we will be able to counteract the averaging 
process to some degree. With the previously generated PE 
filter, we next solved for the data missing region, x m i SS i ng 
that minimizes the following two objective functions, 

Mill f pe X-missing o || 2 + 1 1 a X-missing y || 2 ( 2 ) 

where A is an averaging operator that takes 9 looks, and 
y is an SRTM DEM that covers the missing regions of 
the TOPSAR DEM. Here the x m * S3in g actually has the 
dimensions of the TOPSAR DEM, while the y has the 
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are shown in Fig. 2a and Fig. 2b respectively. The p\ and p 2 
are weights for the two objective functions. 

III. Results 

The cross-validation technique [2] was used to determine the 
optimal weights (Fig. 3) of the two objective functions in (2). 
Since the interpolated region tends to have the same spectral 
content, it is crucial to select a proper window that will be used 
to generate the PE filter. Unless the selected window represents 
the true spectral content of the regions with data missing , it 
is not a good idea to rely heavily on the PE filter. Thus, we 
tried to fit the SRTM DEM as much as possible, allowing 
minimal effect from the PE filter. This can be achieved by 
implementing cross-validation only for the second objective 
function in (2). 

In the case of the example shown in Fig. 2, the minimum 
cross-validation sum of squares (CVSS) was obtained when 
Pi =0.16 and M 2 = 1- Even though only the second objective 
function was used to calculate the CVSS, the mimimum did 
not occur at /ii = 0* because the averaging operator alone 
could not predict the missing data in an adjacent pixel. Fig. 
4 shows the interpolation results applied on Fig. 2, with 
the optimal weighting (b) and two extreme cases (a,c). The 
same method was applied to the entire crater of Sierra Negra 






Fig. 3. Cross-validation sum of squares, ( 2 = 1) 


volcano at Galapagos Islands (Fig. 5). 

IV. Conclusion 

The roll artifact of the aircraft was found to dominate 
in the difference image between the TOPSAR and SRTM 
topographic data This error can be eliminated by subtracting 
the difference from the TOPSAR DEM. Solving the inverse 
problem constrained with a PE Iter and SRTM DEM 
information provided high-quality interpolation results. The 
cross-validation seems to be a proper choice for nding 
the optimal weighting in the inversion. This makes the 
interpolation less biased and guarantees the best ttingto the 
SRTM DEM. The quality of many other TOPSAR DEMs can 
be improved in a similar way. 
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Fig. 4. Interpolated TOPSAR DEM with weights combination of (a) Mi = 1, M2 = 0, (b) Mi = 0.16, M 2 = 1, and (c) Mi = 0, M 2 = 1. Cross-sections 
along A-A’ in (a) are shown in (d). 



Fig. 5. The original TOPSAR DEM (a) and the reconstructed DEM of the same area after interpolation with PE filter and SRTM DEM contraints (b). The 
grayscale is in meters, and the scale is about 12 km across the image. 







